Convective Instability of a Relativistic Ejecta Decelerated by a 
Surrounding Medium: An Origin of Magnetic Fields in GRBs? 

Amir Levinson^ 
ABSTRACT 

Global linear stability analysis of a self-similar solution describing a relativis- 
tic shell decelerated by an ambient medium is performed. The system is shown 
to be subject to the convective Rayleigh- Taylor instability, with a rapid growth 
of eigenmodes having angular scale much smaller than the causality scale. The 
growth rate appears to be largest at the interface separating the shocked ejecta 
and shocked ambient gas. The disturbances produced at the contact interface 
propagate in the shocked media and cause nonlinear oscillations of the forward 
and reverse shock fronts. It is speculated that such oscillations may affect the 
emission from the shocked ejecta in the early afterglow phase of GRBs, and may 
be the origin of the magnetic field in the shocked circum-burst medium. 
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namics, relativity, turbulence 



1. Introduction 

Following shock breakout and a rapid acceleration phase, the relativistic ejecta accu- 
mulated in the explosion starts decelerating, owing to its interaction with the surrounding 
medium. At early times a double shock structure forms, consisting of a forward shock that 
propagates in the ambient medium, a reverse shock crossing the ejecta and a contact inter- 
face separating the shocked ejecta and the shocked ambient medium. In the fireball scenario 
commonly adopted, the naive expectation has been that the crossing of the reverse shock 
should produce an observable optical flash. Despite considerable observational efforts, such 
flashes have not been detected yet. The afterglow emission, which is produced behind the 
forward shock, seem to indicate strong amplification of magnetic fields in the post shock 
region, by some yet unknown mechanism. Moreover, the lightcurve of the afterglow emission 
deviates at early time from that predicted by a simple blast wave model. 
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A question of considerable interest is the stability of the double shock system. Hydrody- 
namic instabilities may lead to strong distortions of the system that may generate turbulence, 
amplify magnetic fields, and affect the emission processes in the afterglow phase. Such effects 
have been studied in the non-relativistic case in connection with young supernovae remnants 
(SNRs). In fact, the idea that the Rayleigh- Taylor (R-T) instability may play an important 
role in the deceleration of a non-relativistic ejecta dates back to Gull (1973), who performed 
ID simulations of young SNRs that incorporate a simple model of convection. Chevalier et 
al. (1992) performed a global hnear stability analysis of a self-similar solution describing the 
interaction of a nonrelativistic ejecta with an ambient medium and found that it is subject to 
a convective instability. They analyzed self-similar perturbations and showed that the flow 
is unstable for modes having angular scale smaller than some critical value. The convective 
growth rate was found to be largest at the contact discontinuity surface and to increase with 
increasing / number of the eigenmodcs. They also performed 2D hydrodynamical simula- 
tions that verified the linear results and enabled them to study the nonlinear evolution of 
the instability. The simulation exhibits rapid growth of fingers from the contact interface 
that saturates, in the nonlinear state, by the Kelvin-Helmholtz instability. Strong distortion 
of the contact and the reverse shock was observed with httle effect on the forward shock. 
Jun & Norman (1996) performed 2 and 3D MHD simulations of the instabihty to study the 
evolution of magnetic fields in the convection zone. They confirmed the rapid growth of 
small scale structure reported by Chevalier et al. (1992), and in addition found strong am- 
plification of ambient magnetic fields in the turbulent fiow around R-T fingers. On average, 
the magnetic field energy density reaches about 0.5% of the energy density of the turbulence, 
but it could well be that the magnetic field amplification was limited by numerical resolution 
in their simulations. The simulations of Chevalier et al. (1992) and Jun Sz Norman (1996) 
support earlier ideas, that the clumpy shell structure observed in young (pre-Sedov stage) 
SNRs such as Tycho, Kepler and Cas A is due to the R-T and K-H instabihty. 

In this paper we extend the linear stability analysis of Chevalier et al. (1992) into the 
relativistic regime. We find that denser ejecta sweeping a lighter ambient gas are subject 
to the R-T instability also in the relativistic case. The stabihty of a double-shock system 
has been investigated by Xiaohu et al. (2002) using the thin shell approximation. However, 
this study is limited to large scale modes and neglects pressure gradients and, therefore, 
excludes the convective instability. Gruzinov (2000) performed a linear stability analysis 
of a Blandford-McKee (BMK) blast wave solution, and found that the BMK solution is 
stable but non- universal, in the sense that some modes decay very slowly as the system 
evolves. Furthermore, the onset of oscillations of an eigenmode of order / has been seen in 
the simulation once the Lorentz factor evolved to F < Z. The conclusion drawn based on 
Gruzinov's findings is that distortion of the shock front at early times may cause significant 
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oscillations during a large portion of its evolution. If the amplitude of these oscillations 
is sufficiently large, and if the same behavior holds in the nonlinear regime then this can 
lead to generation of vorticity in the post shock region (Milosavljevic et al. 2007; Goodman 
& MacFadyen 2007), and the consequent amplification of magnetic fields, as demonstrated 
recently by Zhang et al. (2009). 



2. Analysis 

We consider the interaction of a cold unmagnetized shell with a cold ambient medium 
having a density profile pi = hr~^. The structure formed at early stages consists of a forward 
shock propagating in the ambient medium, a reverse shock propagating in the ejecta and a 
contact discontinuity separating the shocked ambient medium and the shocked ejecta. The 
equations governing the dynamics of the fiow on each side of the contact interface can be 
written in the form, 

, r,(iln7 .,dP dP 

^ln(P/p^)=0, (2) 

P7^(/^7Vt) + VtP = 0, (3) 
at 

where j = is the Lorentz factor of the fiuid, v^- is the tangential component of the 3- 
velocity, which we express as v = f^f + v^-, p, P and h are the proper density, pressure and 
specific enthalpy, 7 is the adiabatic index, and 

r 06 r sm u 0(p 

Equations ([I])-© admit self-similar solutions (Nakamura & Shigeyama, 2006) in cases of a 
freely expanding ejecta characterized by a velocity Ve = r/t at time t after the explosion, 
and a proper density profile 



where 7e = l/yl — is the corresponding Lorentz factor of the unshocked ejecta. The 
Lorentz factors of the forward shock, reverse shock and the contact discontinuity surface, 
denoted by ri(t), r2(t) and Tdt), respectively, evolve with time as Tl = At~"^, Tf = Bt~"^, 
= Ct~"^, with the constants A,B,C and m determined by matching the solutions in 
the shocked ejecta and in the shocked ambient medium at the contact discontinuity. In 
particular m = (6 — 2k) /{n + 2) (Nakamura & Shigeyama, 2006). The velocity of the 



-4- 



ejecta just upstream the reverse shock is Ve{r = = r2{t)/t = 1 — l/[2(m + l)r2], where 
^2(^) = / (1 ~ l/2r2)(it is trajectory of the reverse shock, implying 7g = (m + l)r2. Thus, 
the self-similar solution is applicable only to situations where the ejecta is sufficiently dense, 
such that the reverse shock is non or at most mildly relativistic. 

We have carried out a global linear stability analysis of the self-similar solution outlined 
above. The details will be presented elsewhere (Levinson, in preparation). A preliminary 
account of the method and results is presented below. 

There are total of eight independent variables, four on each side of the contact: P, p, 7, v^. 
The perturbations of these variables were expanded in spherical harmonics. To be more con- 
crete, a perturbed quantity Q {Q = P,p, etc.) is expressed as Q{t,x,0,(f>) = Qo{'tyX)[^ + 
^q(x, t)Yim{0, (p)], where x = {1 + 2(m + l)r^}(l — r/t) is the self-similarity parameter, and 
Qq denotes the unperturbed value. The linearized equations on each side of the contact 
discontinuity were then obtained upon substitution of the perturbed quantities into Eqs. 
dl])-©. Perturbations of the shock fronts and the contact discontinuity of the form 

<5r,(t,^,0) = ^11^(^,0), (6) 
3 

where j = 1,2, c refers to the forward shock, reverse shock, and the contact discontinuity, 
respectively, were assumed. The perturbed shock normals are then n^^ = n^^ + Sn^f^ [k = 
1, 2), with = (— FfcVfc, r^, 0) being the unperturbed normal and 

drik,, = {-Tl6Vk, TlVM, -TkVTSn). (7) 

Here Vk denotes the 3- velocity of the unperturbed shock front and 6Vk = dSr/^/dt. The 
linearized jump conditions at the forward and reverse shocks are given in terms of the 4- 
velocity and the energy-momentum tensor T'^'^ as: 

+ A,(pm'^X] = 0, (8) 
[TrSnj,, + AfcT^-nL] = 0, (9) 

where AkQ = 6Q + {dQo/dr)6rk denotes the Lagrange perturbation of the quantity Q 
at the perturbed surface k. The relations ([8]), ([9]) provide 6 boundary conditions for the 
perturbation equations, 3 on each side of the contact. Two additional boundary conditions 
are obtained from pressure balance and the no flow condition, viz., v — dvc/dt = 0, at the 
contact discontinuity. 

Unlike in the non-relativistic case, the boundary conditions at the shock fronts break 
self-similarity of the perturbations. Specifically, it can be shown that at the forward shock 
surface the perturbation of the tangential velocity is related to the perturbations of the radial 
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velocity and pressure through — 3(^p — ^K)/[2(m + 2)rf] for I ^ 0, with a similar relation, 
though somewhat more involved, at the reverse shock front. Thus, numerical simulations 
of the perturbation equations is needed. The only exception is the spherical mode (/ = 0) 
for which a self-similar solution was obtained analytically. For this solution bQ oc with 
s < 0. It can be shown that one eigenmode of order / = is associated with linear time 
translation of the self-similar solution. For this mode s = —{m + 1). We found another 
eigenmode of order I — that decay somewhat slower. The analytic solution for the I — 
mode has been used both to test the code and as initial condition for the evolution of the 
higher order eigenmodes. Numerical integration of the perturbation equations was performed 
after transforming to the so called Riemann invariants. We identified three variables that 
propagate from the forward shock inwards, three that propagate from the reverse shock 
outwards and two that propagate from the contact, one inwards and one outwards, and 
have chosen the boundary conditions for the Riemann invariants accordingly. It is generally 
found that eigenmodes having an angular scale larger than the horizon scale, specifically 
/(/ + 1) < r^, are stable. Higher order modes are found to be unstable with a growth rate 
that increases with increasing I. An example is exhibited in Fig. 1. The onset of oscillations 
followed by a rapid growth of the initial perturbations is clearly seen. The distortion of the 
contact discontinuity surface becomes nonlinear very early on. The growth is algebraic in 
time t with a growth rate of about 10 in the example shown in Fig. 1; that is, SQ oc t^'^ 
for Q — P, p,v. As seen, the reverse shock responds quickly to the distortion of the contact. 
The forward shock, on the other hand, responds much later, at time when the instability 
near the contact already reached a nonlinear state. 

3. Discussion 

The stability analysis described above seem to indicate strong convective instability 
at early stages of the evolution of a dense cjccta as it sweeps a lighter ambient gas. The 
growth rate appears to be largest at the contact discontinuity and for higher order modes. 
Disturbances at the interface separating the shocked cjccta and the shocked ambient medium 
propagate away from the contact discontinuity and cause nonlinear distortions of the shock 
fronts. The reverse shock responds quickly to the distortion of the contact. Propagation 
of the signal to the forward shock is much slower. At any rate, the instability near the 
contact becomes nonhnear well before the signal arrives at the forward shock, so full MHD 
simulations are needed to resolve the effect of the instabihty on the forward shock. It is 
naively expected that the instability will be strongly suppressed in cases where the ejecta is 
highly magnetized and/or if the reverse shock is highly relativistic. On the other hand, if 
the magnetic field strength in the unshocked ejecta is smaller than that required to suppress 
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the instability but still much larger than that of the ambient medium, then at early stages 
mixing of the magnetized ejecta with the shocked ambient gas via growth of R-T fingers 
can give rise to a strong amplification of the magnetic fields behind the forward shock. Full 
simulations are required to quantify the conditions under which the instability is effective. 

The nonhnear distortions of the contact and the shock fronts should generate turbulence 
in the shocked fluids on both sides of the contact discontinuity. At early stages this may 
strongly affect particle acceleration and the emission processes. It is tempting to speculate 
that the lack of observed optical flashes, that are anticipated in the "standard" model, and 
the fact that the early afterglow emission observed in many sources is inconsistent with the 
prediction of the blast wave model may be attributed to the instability discussed here. In 
any case, it is clear that a careful analysis that takes account of this process is required to 
better understand the observational characteristics of the emission during the early afterglow 
phase. 

The stability analysis of the Blandford-McKee solution performed by Gruzinov (2000) 
suggests that it may be a very slow attractor. Linear perturbations of the forward shock 
in the B-M phase decay very slowly. Whether this behavior continues also in the nonlinear 
regime is unclear yet. If it does then it is anticipated that the growth of R-T fingers and, 
perhaps, nonlinear oscillations of the forward shock itself that are induced by the convective 
instability may be a source of vorticity during a long portion of the evolution of the blast 
wave. As demonstrated recently by Zhang ct al. (2009) the induced turbulence can amplify 
weak magnetic fields. Their simulation seem to converge at a saturation level of ~ 
5 X 10^'^, weakly dependent on the initial magnetic field strength. This process may provide 
an explanation for the origin of the strong magnetic fields inferred behind the coUisionless 
shock in the afterglow phase. 

Unfortunately, the linear analysis outlined above is restricted to a limited set of condi- 
tions under which the unperturbed self-similar solution of Nakamura & Shigeyama (2006) is 
applicable. Full 3D MHD simulations are required to study this process in other situations, 
and to follow the evolution of the convective instability in the nonlinear state. As illustrated 
above, high resolution simulations that can resolve angular scales << l/F are required, 
posing a great numerical challenge. We believe that our findings strongly motivate such 
efforts. 

I thank A. Ditkowski for a technical help in the development of the code, and M. 
Alloy, A. MacFadyen, E. Nakar and E. Waxman for enlightening discussions. This work was 
supported by an ISF grant for the Israeli Center for High Energy Astrophysics, and by the 
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Fig. 1. — Time evolution of the perturbations for n = 1.1, k = 2 and /(/ + l)/r^(to) = 10^, 
here ri(to) is the initial Lorentz factor of the forward shock. Upper panels: perturba- 
tion of the contact discontinuity surface (left panel) and relative pressure perturbation, 
5P{t)/6P{to), of the shocked ambient medium at the contact (right panel). Bottom panels: 
perturbations of the reverse shock surface (left), and relative pressure perturbation at the 
reverse shock front (right). The initial perturbations of the contact discontinuity and the 
reverse shock surface in this example are Srjrc = 6 x 10~^ and Srs/vg — 10~^, respectively. 



